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I.  Introduction 


For  over  thirty  years  Mathematical  Programming  Systems  (MPS’s)  have  been  built  arround  sophisticated 
implementations  of  the  simplex  method,  with  its  many  variations,  refinements  and  extensions.  During 
almost  this  entire  period  the  simplex  method  was  unchallenged  as  the  method  of  choice  for  solving  linear 
programs,  and  by  extension  those  other  classes  of  problems,  such  as  mixed  integer  and  some  kinds  of  nonlinear 
programs,  which  use  sequences  of  linear  programs  (LP’s).  In  the  light  of  new  developments  in  mathematical 
programming,  this  situation  now  requires  re-examination. 


The  projective  method,  published  by  Karmarkar  (1984),  is  an  interior  point  method,  which  at  first  sight  has 
nothing  in  common  with  the  simplex  method,  proceeding  as  it  does  through  the  (relative)  interior  of  the 
feasible  regi  i,  rather  than  from  vertex  to  vertex  of  the  polytope  defining  that  region.  While  claims  have 
been  made  fc  the  theoretical  virtues  of  other  methods  (notably  the  ellipsoid  method  of  Kachiyan(1979)), 
Karmarkar’s  v  rk  has  drawn  much  attention  because  of  widely  publicised  claims  that  his  method  is  orders 
of  magnitude  jc.’ter  than  the  simplex  method  and  some  of  the  commercial  MPS’s  that  employ  it. 


Much  of  the  controversy  remains  to  be  resolved.  It  turns  out  that,  viewed  in  the  geometry  in  which  Dantzig 
originally  concievrd  the  simplex  method  (see  Dantzig(l963)),  the  two  methods  are  not  wholly  unrelated 
(Stone  and  Tovey(1986)).  Even  more  intriguingly,  it  has  been  shown  by  Gill  et  al  (1985)  that  the  projective 
method  is  equivalent  to  a  special  case  of  a  barrier  function  method  (see  Fiacco  and  McCormick(1968)). 
This  has  lead  to  a  much  better  understanding  of  the  mathematical  and  computational  proceedures  involved, 
which  is  now  being  backed  up  by  independant  computational  results.  In  particular,  Gill  et  al  (1985)  obtained 
computational  results  showing  the  barrier  method  to  be  comparable  in  efficiency  with  the  simplex  method. 
Subequently,  Adler  et  al  (1986)  presented  computational  results  which  indicate  that  an  affine  form  of  the 
projective  method  (see  Vanderbei  et  al  (1985)),  operating  on  the  dual  LP  can  be  faster  than  some  simplex 
implementations  by  a  factor  of  two  or  more  on  many  problems. 


f 


Even  though  these  computational  results  are  very  preliminary,  they  are  impressive  enough  to  warrant  a 
realistic  evaluation  of  interior  point  methodologies  in  the  context  of  the  large  scale  MPS’s  actually  used  to  — 
solve  large  scale  real  world  problems.  Such  integration  and  testingjs  the  topic  of  this  paper. 


□  □ 


2.  Special  Purpose  Codes  and  MPS’s 


Development  of  nev  mathematical  programming  algorithms  (usually  for  problems  of  particular  structure) 
has  frequently  lead  to  the  development  of  special  pupose  codes,  usually  written  in  FORTRAN  or  some  other 
high  level  language.  These  are  often  used  to  test  the  algorithms  and  then  *  cleaned  up”  for  use  on  practical 
problems.  However,  in  most  cases  the  facilities  provided  are  rudimentary  •  the  ability  to  input  a  model  in 
standard  form  (such  as  MPS  format),  an  optimizer,  and  solution  output  (perhaps  also  in  NIPS  format)  for 
the  best  solution  obtained.  It  is  rare  to  find  facilities  for  revising  and  restarting  models,  parametrics  or 
ranging,  or  any  of  the  other  "bells  and  whistles*  provided  by  a  large  scale  MPS  (see  Orchard-Hays  (1968) 
for  a  description  of  MPS  facilities  and  terminology).  There  are,  of  course,  exceptions  which  prove  the  rule; 
notably  the  restart  and  model  revision  tools  provided  in  MINOS  (see  Murtagh  and  Saunders  (1985))  and 
the  "X.  System"  of  Brown  and  Graves  (1983). 

A  large  scale  MPS  is  a  complicated  software  system  which  may  be  dominated  not  by  algorithmic  considera¬ 
tions,  but  by  model  management,  sensitivity  analysis,  case  study  and  reporting  requirements.  The  end  user 
usually  does  not  know  (or  care)  how  the  results  were  produced,  but  often  requires  voluminous  reports  and 
the  ability  to  work  from  a  large  database.  The  resulting  environment  may  then  be  something  like  that  in 
Figure  1;  a  simplified  schematic  of  Ketron  Management  Science’s  proprietary  MPSHI  system. 

The  data  path  on  the  right  of  Figure  1  uses  time-honored  MPS  card-image  input,  which  is  CONVERTed  to 
an  intermediate  packed  form  on  the  PROBFUE.  This  is  followed  by  SETUP,  which  selects  appropriate  right 
hand  side,  objective  and  bound  information  for  the  specified  case  and  creates  a  packed "workfile"  or  "SETUP 
problem".  In  MPSHI  (as  in  all  descendants  of  MPS/360)  this  is  an  out-of-core  file.  (See  also  Benichou  et 
*1(1977)). 

The  data  path  on  the  left  of  Figure  1  uses  the  DATAFORM  model  management  language  to  hold  the  database, 
build  models,  save  solution  cases,  create  reports,  perform  case  studies,  etc.  Here  models  are  saved  in  tree- 
structured  form  on  the  random-access  ACTFILE,  which  contains  all  data  pertinant  to  them,  rather  than 
on  a  PROBFUE  (though  conversion  is  possible  between  the  two  forms).  The  Program  Control  Language 
(PCL)  verb  READY  performs  analagously  with  SETUP  in  producing  a  workfile  from  the  ACTFILE. 

The  workfile  produced  by  these  data  processing  steps  is  the  model  representation  to  be  operated  on  by  the 
algorithmic  modules.  These  include  the  PRIMAL  or  VARIFORM  algorithm  (with  the  option  of  using  GUB), 
DUAL  and  numerous  Ranging  and  Parametric  algorithms.  Solution  reporting,  some  forms  of  matrix  editing, 
reoptimization  from  an  old  basis  and  many  other  functions  must  also  be  done  at  the  workfile  leveL 

MPSHI  (and  some  other  systems)  may  carry  the  above  process  a  step  further  and  repack  the  workfile 
to  give  an  even  more  compact  representation  using  the  concept  of  "super-sparsity"  (see  Kaian(1971)  or 
Greenberg(1978)).  This  concept  makes  use  of  the  fact  that  the  number  of  "unique  values"  even  in  a  sparse 
matrix  may  be  relatively  small,  and  in  particular,  many  of  them  may  be  +1  or  -1.  This  may  be  taken 
advantage  of  by  using  a  "pool"  of  unique  values,  which  is  referenced  by  the  individual  nonzeros  and  by 
avoiding  multiplications  by  unit  coefficients.  Thus  the  standard  workfile  represents  vectors  (columns)  as  in 
Figure  2(a),  while  WHIZARD,  our  high-speed  in-core  optimizer  uses  a  representation  as  Figure  2(b).  If 
WHIZARD  is  used  for  rapid  optimization,  the  solution  information  (which  variables  are  basic  and  which 
at  bound)  must  be  returned  to  the  model  in  standard  workfile  format,  where  we  can  then  INVERT  and 


recompute  the  solution  in  the  required  form  for  use  by  other  modules. 


It  should  be  clear  that  such  MPS’s  posses  a  great  deal  of  power  and  flexibility  which  their  developers  and 
users  are  loath  to  lose  or  be  forced  to  recreate.  The  question  of  when  and  if  some  algorithm  other  than  the 
revised  simplex  method  (around  wh'ch  such  systems  were  designed)  should  be  incorporated  is  thus  a  serious 
one.  Beale  and  Tomlin  (1970)  elucidated  three  criteria  for  incorporating  new  algorithms: 


(a)  They  are  reasonably  easy  to  implement. 

(b)  They  are  reasonably  easy  to  use. 

(c)  They  enable  a  significant  class  of  problems  to  be  solved  more  easily  in  an  MPS  framework. 

Good  examples  satisfying  these  criteria  have  been  GUB  (Generalized  Upper  Bounding)  and  Branch  and 
Bound  methods.  Given  a  super-sparse  capability  of  the  WHIZARD  type,  Tomlin  and  Welch  (1984, 1985) 
showed  that  special  primal  simplex  algorithms  for  pure  and  generalized  networks  were  also  excellent  examples. 
On  the  other  hand,  attempts  to  incorporate  decomposition  as  an  integral  part  of  MPS’s  (as  opposed  to  user 
implementation  -  see  Ho  and  Loute  (1983))  have  been  less  successful. 

When  we  come  to  examine  the  newer  interior  point  methods  in  this  light,  we  find  that  criterion  (b)  is  easily 
dealt  with  *  the  new  algorithm  is  simply  called  in  place  of  the  old  for  primal  optimization.  Criterion  (a)  is 
the  topic  of  most  of  the  rest  of  this  paper.  Criterion  (c)  is  of  course  still  a  matter  of  major  controversy,  but 
we  have  preliminary  computational  results  which  indicate  positive  results  for  at  least  one  significant  class  of 
problems. 


3.  MPS  Data  Structures 


We  have  already  described  some  of  the  data  structures  used  for  matrix  representation  at  the  algorithmic 
level  in  MPSIII.  We  should  also  consider  the  related  data  structures  used  in  the  real  work  of  the  revised 
simplex  method.  These  are  the  structures  of  the  "transformations"  (sometimes  known  as  "etas”)  used  to 
update  vectors.  The  basis  (or  its  inverse)  is  factorized  into  a  product  of  elementary  transformations: 

B-x^Ek....E2Ei  (3.1) 


and  is  used  to  compute: 

{flit  •••»  =  E  (fly,  1 •••)  o/»}  (3.2) 


ft 


Note  that  usually  several  vectors  are  transformed  at  once  (see  Orchard-Hays  (1968)  and  Forrest  and  Tomlin 
(1972)).  The  nonunit  column  (or  row)  of  the  elementary  transformations  may  be  stored  analagously  with 
the  way  in  which  matrix  columns  are  stored  for  the  sparse  and  super-sparse  structures. 

Finally  we  indicate  the  way  in  which  memory  is  allocated  in  the  standard  MPSIH  SETUP  mode  and  in 
WHIZARD  in  Figure  3.  The  WHIZARD  structure  is  inherantiy  more  flexible,  and  is  the  one  we  have  used 
for  other  special  algorithm  implementations.  Note  that  we  have  shown  a  larger  partition  for  the  WHIZARD 
structure,  since  the  WHIZ1N  conversion  module  seeks  out  as  large  a  work  space  as  possible.  The  size, 
manageability  and  availability  of  this  space  (normally  used  for  transformations)  is  important  for  interior 
point  methods  since  they  use  a  number  of  arrays  of  the  column  dimension  of  the  LP.  Allocation  of  such 
arrays  in  the  SETUP  mode  would  be  very  difficult. 

4.  Choice  of  Method(s) 

Since  the  original  appearance  of  Karmarkar’s  work  there  have  been  numerous  variations  and  extensions,  and 
a  choice  (or  choices)  for  implementation  must  be  made.  We  chose  the  Newton  barrier  method  (see  Gill  et  al 
(1986))  for  our  initial  implementation  for  the  following  reasons: 

(a)  Familiarity,  due  to  our  collaboration  on  the  cited  reference. 

(b)  The  fact  that  Karmarkar’s  projective  method  and  the  affine  method  of  Vanderbei  et  al  (1985)  may  be 
extracted  as  special  cases  of  this  method. 

(c)  The  ability  to  handle  problems  in  standard  form,  with  bounds  and  ranges,  with  reasonable  efficiency. 

(d)  The  fact  that  it  works  with  a  primal  feasible  solution  (if  there  is  one).  This  allows  for  early  feasible 
termination. 

We  shall  not  recapitulate  the  details  of  the  algorithm  (see  GUI  et  al  (1986)),  merely  give  an  outline  which 
emphasizes  the  critical  sub-algorithms  so  that  we  see  how  they  may  affect  the  MPS.  The  problem  is  stated 
as: 

wn'jT{c}z,  - nhz}) 
i 

subject  to : 

Az  =  b 
0  <Zj<Uj 

where  p  — ♦  0. 

To  outline  the  algorithm  let  us  define: 

D  = 


Then  the  steps  of  an  iteration  are: 


(1)  If  p  and  ||r||  are  sufficiently  small  -  STOP. 

(2)  If  "appropriate"  reduce  p  and  recompute  r. 

(3)  Solve  a  least  squares  problem: 

min  ||r  -  DAr$*j|.  (4.1) 

(4)  Update  the  "pi  values*  and  "reduced  costs": 

f  +  i^-d-AT6t} 
and  compute  the  search  direction  p  as: 

r  =  Dd-  pe,  p  =  -(l/p)Dr. 

(5)  Calculate  the  steplcngth  a. 

(6)  Update  x  *-  x  +  ap.  GO  TO  (1). 

The  most  critical  step  in  each  iteration  is  the  gradient  projection,  that  is  solving  the  least  squares  problem 
(4.1).  As  in  Gill  et  al.(l986)  we  use  a  preconditioned  conjugate  gradient  method  (the  LSQR  method  of  Paige 
and  Saunders  (1982)).  This  involves  the  Cholesky  factorization  of  (an  often  approximate)  normal  equation 
matrix 

PLLTPT  =  AD2Xr:=AD2AT, 

where  the  permutation  matrix  P  is  chosen  so  that  the  lower  triangular  factor  L  is  sparse.  We  then  solve  the 
preconditioned  problem 

mm  ]jr  -  DATPL~Ty\\  (4.2) 

and  recover  St  -  PL~Ty. 

The  essential  steps  in  LSQR  are  calculations  of  the  form 

a  *-  u  +  DAT{PL~Tv))  (4.3) 

and 

tM-  v  +  L~lPT(ADu).  (4.4) 


5.  Some  Implementation  Considerations 

The  WHIZARD  environment  turns  out  to  be  quite  convenient  for  implementation  of  the  Newton  barrier 
method  in  the  (assembly  language)  MPSIII  system.  A  first  advantage  is  that  all  of  the  Presolve/Postsolve 
facilities  of  WHIZARD  are  automatically  available  to  reduce  and  simplify  the  model,  in  the  same  way  as 
they  are  for  WHIZNET  (see  Tomlin  and  Welch  (1983,1985)).  The  existing  WHIZIN  module  also  handles 
memory  management  and  allocation  effectively.  In  particular  most  of  the  arrays  automatically  assigned 
when  WHIZIN  reconfigures  memory  for  the  simplex  method  have  immediate  analogues  or  obvious  uses  in 
the  barrier  context: 

(1)  The  t  vector  is  used  in  both  algorithms. 

(2)  The  p  region  usually  used  for  the  values  of  the  basic  variables  holds  the  right  hand  side. 

(3)  The  "basis  headings”  array  or  "H-region"  is  available  for  saying  the  row  permutation  P  and  its  inverse. 

(4)  The  7  row-length  work  regions  are  available  to  hold  the  "artificial  column"  for  phase  I,  h,  and  the  work 
vectors  for  LSQR. 

The  large  transformation  block  may  be  segmented  to  store: 

(5)  The  column-size  vectors  z,c,d,r. 

(6)  The  preconditioner  L. 

There  is  one  more  large  data  structure,  which  allows  access  to  A  row-wise  so  that  AD2AT  may  be  computed. 
This  is  explained  below. 

An  important  decision  had  to  be  made  on  the  extent  to  which  we  would  choose  data  structures  which 
accomodate  existing  subroutines  and  make  implementation  easier,  particularly  in  the  critical  steps  of  com¬ 
puting  and  applying  the  preconditioner.  A  great  deal  of  the  simplex  code  mid  be  "mined"  for  the  barrier 
method.  For  example,  note  that  the  computational  steps  in  (4.3)  consist  of  updating  t;  by  L~T ,  equivalent 
to  a  BTRAN  of  a  pricing  vector  in  the  simplex  method  (3.3),  and  multiplication  of  DAT  by  the  resulting 
vector  involves  essentially  the  same  work  as  a  full  (scaled)  PRICE  of  the  matrix.  Similarly  in  (4.4)  we  see 
that  forming  ADu  requires  essentially  the  same  work  as  performing  a  CHECK  for  a  (scaled)  solution  u. 
The  update  of  the  permuted  ADu  by  L“l  corresponds  to  the  FTRAN  of  a  vector  in  the  simplex  method 
(3.2).  The  computation  of  the  Cholesky  factors  LlTthemselves  could  be  regarded  as  a  special  symmetric 
application  of  the  sparse  factorization  used  in  simplex  INVERT  routines. 

It  is  a  natural  temptation  to  use  the  WHIZARD  transformation  data  structures  and  modify  INVERT, 
FTRAN,  BTRAN  etc.  to  exploit  existing  code.  We  decided  not  to  do  this  for  a  number  of  reasons: 

(a)  W<  expect  L  to  be  not  at  all  super-sparse,  unlike  the  basis  factors  produced  by  INVERT  for  the  simplex 


method,  and  indeed  to  be  much  denser  altogether. 

(b)  The  simplex  FTRAN  and  BTRAN  (see  (3.2-3))  routines  are  designed  with  the  efficient  update  of  multiple 
vectors  in  mind.  We  do  not  wish  to  do  that  at  this  stage. 

(c)  Published  algorithms  are  available  (see  George  and  Liu  (1981))  which  can  be  easily  integrated,  using 
simple  data  structures,  and  which  predict  storage  requirements  for  L,  thus  enabling  us  to  avoid  exceeding 
space  available. 

We  thus  chose  to  adopt  the  data  structures  used  by  George  and  Liu  (1981),  including  the  use  of  the  "com¬ 
pressed  storage  scheme",  with  the  modification  that  the  row  index  array  is  of  full  word  length  so  that  the 
indices  multiplied  by  eight  may  be  stored  (to  avoid  shifts  in  inner  loops).  Pointers  into  the  index  and  value 
arrays  are  also  shifted  appropriately. 

6.  Auxilliary  Processing 

There  are  several  non-trivial  proceedures  which  must  be  carried  out  before  the  algorithm  can  be  even  begun. 
In  particular,  the  permutation  P,  the  nonzero  structure  of  L  and  a  representation  of  AT  must  be  determined. 

To  determine  P  we  require  the  non-zero  structure  of  AAT.  We  do  this  by  processing  the  columns  of  A 
sequentially  and  building  a  list  of  index  pain  (11,1*2)  in  each  column,  each  pair  in  a  full  word.  When  all 
columns  have  been  processed,  or  periodically  if  array  space  is  exceeded,  this  list  b  sorted  (using  Shell  sort). 
Duplicate  pairs  may  then  be  purged  in  a  single  pass,  and  the  data  reconfigured  into  a  form  suitable  for  a 
minimum  degree  ordering,  which  gives  P  (see  Liu  (1985)).  This  'is  then  followed  by  a  symbolic  factorization, 
derived  from  George  and  Liu  (1981).  As  stated  above,  the  data  structure  is  modified  for  efficient  application 
of  the  preconditioner  L.  Once  the  nonzero  structure  of  L  is  available,  the  memory  requirements  for  the 
algorithm  can  be  computed,  and  if  insufficient  space  is  available  the  proceedure  may  be  halted. 

The  remaining  major  data  structure  is  required  so  that  we  may  compute 

LLT  =  PTAD2ATP  (6.1) 

where  A  and  D  are  modified  to  exclude  dense  columns  and  those  corresponding  to  tiny  r;.  Thb  requires 
either  a  much  more  expensive  updating  proceedure  for  the  data  structure  L,  if  we  access  A  only  by  column, 
or  the  ability  to  access  the  matrix  row-wise.  We  chose  the  latter  alternative.  The  data  structure  used  to  do 
thb  is  as  follows: 

Since  we  already  have  a  count  of  the  number  of  nonzeros  in  each  row  from  earlier  processing  it  is  easy  to  set 
up  an  array  which  contains  a  single  word  for  each  useful  nonzero  (i.e.  those  not  eliminated  by  PRESOLVE) 
in  the  matrix.  A  row  length  array  points  to  the  first  entry  for  each  row  (see  Figure  4)  in  permutation 
order.  Each  entry  in  the  ROWLNK  array  consbts  of  a  3-byte  column  number  for  the  element  and  a  1-byte 
offset  which  gives  the  offset  of  the  nonzero  in  the  super-sparse  packed  WHIZARD  matrix.  There  b  also 
a  column-length  COLADR  array  (also  used  for  other  purposes)  which  gives  the  address  of  each  column  in 
the  matrix  and  a  1-byte  offset  which  specifies  where  the  non-unit  elements  begin  relative  to  the  top  the 


column  description.  The  offsets  in  the  two  arrays  may  thus  be  compared  logically  to  see  if  a  nonzero  in  the 
row  being  referenced  is  a  unit  entry  or  not,  and  treat  it  accordingly.  Note  that  the  ROWLNK  array,  once 
allocated,  can  be  built  in  a  single  pass  of  the  matrix,  using  only  row-size  work  arrays.  We  should  also  point 
out  that  the  1-byte  offsets  limit  this  structure  to  columns  represented  by  at  most  256  bytes.  However,  this 
is  no  restriction,  since  such  columns  would  contain  over  120  nonzeros,  and  we  do  not  wish  to  consider  such 
columns  when  computing  the  preconditioner  anyway. 

7.  Transition  to  a  Basic  Solution 

No  matter  how  rapidly  an  interior  point  algorithm  may  perform,  it  will  not  usually  lead  to  a  in  instantly 
identifiable  optimal  basic  solution.  However  in  many,  if  not  most  applications  a  basic  solution  (or  at  least 
one  with  a  minimal  number  of  nonzero  values)  is  very  important.  Reasons  for  this  are: 

(1)  It  is  desireable  in  practice  to  keep  the  number  of  "active”  activities  small. 

(2)  Basic  solutions  are  more  "nearly  integer”  in  many  models,  e.g.  production-distribution  models. 

(3)  It  is  easier  to  save  and  restore  model  status. 

(4)  All  the  standard  postoptimality  and  sensitivity  analysis  procedures,  including  ranging  and  parametrics, 
assume  a  basic  optimum. 

(5)  Use  of  nonlinear  procedures  requiring  repeated  solution  of  modified  LP  problems  (e.g.  MIP  by  branch 
and  bound  and  SLP)  seem  to  remain  dependant  on  variants  of  the  simplex  method. 


The  idea  of  using  some  non-simplex  method  to  perform  part  of  the  optimization  work  and  then  switching  to 
the  simplex  method  is  by  no  means  new.  In  fact  a  "BASIC"  technique  for  purifying  non-basic  solutions  has 
been  a  part  of  MPS’s  for  many  years  (it  appears  in  the  early  MPS/360  documentation,  see  also  the  MPSIH 
User’s  Manual).  Among  the  uses  that  have  been  made  of  this  proceedure  are  transitions  from  Dantzig-Wolfe 
decomposition  (see  Ho  and  Tomlin  (1977)),  from  an  Augmented  Lagrangian  method  (see  Beale,  Hattersley 
and  James  (1985))  and  projective  /barrier  methods  (see  Gill  et  al.(1986)).  However,  in  the  latter  case  the 
transition  was  carried  out  through  an  external  interface,  which  is  inefficient  and  clumsy.  The  clear  remedy 
is  to  implement  a  BASIC  procedure  internally  in  WHIZARD  to  process  the  nonbasic  solution.  We  realized 
this,  but  made  an  attempt  to  avoid  it. 


Our  experiments  with  interior  point  methods  had  lead  to  the  observation  that  there  were  invariably  less 
than  m  solution  values  significantly  away  from  zero  (or  their  bounds).  In  other  words,  our  problem  is  not  an 
excess  of  candidates  for  basic  status.  It  therefore  seemed  possible  that  a  simple  classification  and  pivoting 
scheme  would  suffice. 


Let  us  assume  some  general  "significance  tolerance"  t  and  suppose  we  have  obtained  "optimal"  primal  and 
dual  values  *  and  r .  The  latter  are  used  to  calculate  the  reduced  costs  (d;-  values)  and  we  classify  the  columns 
and  rows  as  in  Table  1.  The  action  to  be  tried  for  relevant  pairs  of  categories,  based  on  complementarity 
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Cols 


Rows 

- GO 

n, 

PIVOT 

HERE 

SECONDARY 

PIVOT 

FIX 

TO  ZERO 

Slack<  £ 

W<£ 

SECONDARY 

PIVOT 

TERTIARY 

PIVOT 

Slack>  £ 

W  < « 

BASIC 

SLACK 

Table  1 


considerations,  is  given  in  the  body  of  the  table. 

The  primary  attempt  is  to  pivot  columns  with  positive  values  on  active  constraints  (with  x*  ^  0).  Normally 
there  would  be  an  imbalance  of  these,  which  requires  a  secondary  choice  of  either  pivoting  columns  with 
positive  zy  on  rows  with  x<  :=  0,  or  columns  with  a  zero  Zj  on  active  constraints.  There  may  then  be  some 
tertiary  pivots  necessary  to  use  up  the  remaining  rows  which  are  not  assigned  basic  slacks.  If  all  went  well, 
this  scheme  should  lead  to  an  optimum  basic  feasible  solution. 

In  practice  this  naive  scheme  almost  always  failed.  This  seems  to  be  because,  even  for  very  sparse  and 
degenerate  problems,  there  is  significant  linear  dependance  between  the  columns  corresponding  to  positive 
zy  in  the  solutions  produced  by  the  interior  point  method.  Attempts  to  pivot  these  columns  into  a  basis 
therefore  are  doomed  to  failure,  and  we  are  in  the  same  position  as  when  a  singular  ’’basis”  is  provided  by 
a  user  *  columns  must  be  dropped,  and  the  starting  solution  is  infeasible.  The  number  of  simplex  iterations 
needed  to  recover  feasibility  and  then  optimality  may  be  large. 

Given  the  failure  of  the  naive  scheme,  there  was  no  alternative  to  implementing  a  version  of  the  BASIC 
algorithm  within  WHIZARD  itself.  The  steps  of  this  algorithm  are  not  widely  known,  so  we  briefly  restate 
them  here: 

Suppose  all  variables  are  at  tentative  values  zj-  and  that: 

Ax*  =  6, 0  <  zj-  <  Uj 

The  normal  simplex  iteration,  pivoting  on  row  r,  computes: 

ft  <-  ft  -  <»ij ft  (»  /  r) 

when  column  j  enters  the  basis.  .  assume  that,  as  usual,  if  a  vector  has  dy  >  0  it  will  be  decreased  (as  if  it 
was  coming  in  from  it's  upper  bound).  In  this  case  -ay  will  be  FTRANned  and  used  in  CHUZR.  If  dy  <  0 
it  is  to  be  increased  in  the  usual  way.  The  explicit  calculations  arc 


If  fa  <  ( Uj  -  z*)  do  a  regular  pivot,  i.e.  modify  P  by  0rary.  Set  fa- fa  +  x*. 

If  fa  >  ( U y  -  Zj)  do  a  bound  flip  type  move.  Le.  modify  $  by  (Uj  -  z*)cry  and  mark  zy  as  being  at  it’s 
upper  bound  (of  Dy). 

Case(2)  dy  >  0 

If  8r  >  x\  do  a  bound  flip  type  move.  i.e.  modify  0  by  z‘ay  and  mark  zy  as  being  at  it’s  lower  bound 
(of  0)7 

If  8t  <  zj-  do  a  regular  pivot.  i.e.  modify  0  by  0rory.  Set  fa  =  z*j  -  fa. 

Depending  on  the  starting  nonbasic  solution  there  is  no  guarantee  that  the  result  will  be  feasible  or  optimal, 
only  that  the  sum  of  infeasibilities  or  objective  will  not  be  degraded.  It  is  therefore  necessary  to  call  the 
simplex  method  to  verify  or  achieve  optimality  before  entering  Postsolve. 


8.  Preliminary  Computational  Results 

It  is  not  our  intention  to  give  detailed  or  extensive  computational  experience  in  this  paper  (this  will  be 
done  in  a  later  publication).  Some  prelimimy  results  are  of  interest,  however.  Table  2  gives  the  original 
and  reduced  (Presolved)  dimensions  of  a  subset  of  the  problems  in  Gill  et  aJ.(1986).  Note  the  significant 
reductions  achieved  for  some  models,  particularly  the  highly  degenerate  NZFRJ  forestry  model.  If  the  large 
structurally  degenerate  part  of  the  model  is  not  removed  in  this  class  of  models  the  behaviour  of  the  simplex 
method  can  be  badly  affected,  even  using  some  anti-degeneracy  procedure,  whereas  we  see  in  Table  3  that 
the  WHIZARD  "fast  primal"  starting  from  a  "crashed"  basis  easily  outperforms  the  current  interior  point 
implementation.  This  is  true  of  the  other  models  also,  with  one  exception. 


Problem  Name 

Original  Dimensions 

Reduced  Dimensions 

Rows 

Columns 

Rows 

Columns 

AFIRO 

28 

32 

28 

32 

SHARE2B 

97 

79 

97 

79 

BRANDY 

221 

249 

131 

222 

BANDM 

306 

472 

244 

322 

DEGEN2 

445 

534 

442 

534 

NZFRI 

624 

3521 

415 

2034 

Tible  2 


The  model  DEGEN2  is  also  highly  degenerate,  but  does  not  exhibit  a  significant  structurally  degenerate, 
removable  component  of  the  model.  The  simplex  method  requires  a  great  many  iterations  for  this  class  of 
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models  and  the  barrier  method  is  able  to  win  by  a  significant  margin,  even  with  quite  a  lengthy  BASIC 
and  simplex  phase  added  on  to  achieve  a  basic  optimum.  We  have  therefore  found  at  least  one  class  of 
models  for  which  interior  methods  show  promise.  More  generally,  we  should  point  out  that  it  seems  that 
many  of  the  models  used  elsewhere  in  comparisons  to  show  interior  methods  to  advantage  seem  to  be  highly 
degenerate.  Furthermore,  most  comparisons  have  been  made  with  MINOS,  which  is  a  FORTRAN  code  with 
no  anti-degeneracy  procedure  and  no  Presolve  to  remove  structural  degeneracy.  Our  comparisons,  on  the 
other  hand,  are  with  a  state-of-the-art  assembly  language  simplex  code,  with  both  algorithms  in  the  same 
environment. 


Problem  Nzme 

WHIZARD 
FAST  PRIMAL 

WHIZARD 

BARRIER 

BARRIER  PLUS 
BASIC 

AFIRO 

0.06 

0.018 

0.020 

SHARE2B 

0.36 

0.96 

1.14 

BRANDY 

1.32 

3.72 

4.50 

BANDM 

1.86 

5.28 

7.02 

DEGEN2 

35.58 

18.54 

24.00 

NZFRI 

6.06 

20.16 

29.64 

Table  3 


9.  Further  Work 

The  implementation  and  results  we  have  described  represent  only  a  beginning,  with  the  need  for  a  number 
of  improvements  recognized.  The  efficiency  of  the  WHIZARD  internal  BASIC  procedure  can  certainly  be 
improved.  Almost  certainly,  other  improvements  can  be  made  to  the  interior  point  method  itself.  We  have 
not  yet  tried  the  straightforward  Vanderbei  et  al.(1985)  algorithm  or  the  approximate  trajectory  method 
used  in  Adler  et  al.(1986).  Similarly  there  are  possible  improvements  to  be  gained  by  using  different  sparse 
factorization  techniques  along  the  lines  given  by  Gustavson  et  al.(1970).  The  possibility  of  updating  the 
LLT  factors,  as  suggested  by  Karmarkar  (1984)  and  Shanno  (1985)  has  yet  to  be  evaluated,  as  have  dual 
algorithms  (Osborne  (1986)). 

A  pressing  question  for  general  use  of  interior  methods  in  MPS’s  is  whether  restart  procedures  for  modified 
models  can  be  found  which  compete  with  variants  of  the  simplex  method.  Until  this  question  is  answered 
there  is  no  alternative  to  using  simplex-ba3ed  technology  for  such  procedures  as  branch  and  bound,  sequential 
LP,  or  even  regular  production  use  of  many  large-scale  IP  models. 
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10.  Conclusion 


Despite  the  dramatic  differences  in  philosophy  and  approach,  it  seems  that  interior  point  LP  methods  can 
be  implemented  with  reasonable  efficiency,  and  even  elegance,  in  production  Mathematical  Programming 
Systems  based  on  the  simplex  method.  What  is  not  yet  known  is  what  effect  this  will  have  on  the  economical 
solution  of  which  classes  of  large-scale  LP  models.  We  have  tentatively  identified  one  class,  and  certain  other 
models  with  special  structure  (e.g.  block  angular  and  staircase)  seem  good  candidates.  We  expect  to  report 
further  on  implementation  and  computational  experience  in  subsequent  papers. 
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